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Abstract 

A probabilistic method for solving time-dependent load-transfer models of 
fracture is developed. It is applicable to any rule of load redistribution, i.e, 
local, hierarchical, etc. In the new method, the fluctuations are generated 
during the breaking process (annealed randomness) while in the usual method, 
the random lifetimes are fixed at the beginning (quenched disorder). Both 
approaches are equivalent. 
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I. INTRODUCTION 



The modelization of fracture in disordered systems is a subject of great interest in natural 
and artificial materials 0. A time-dependent method to describe the failure of materials 
under stress, within the fiber-bundle paradigm, was proposed by Coleman M . In this model a 
set (bundle) of elements (fibers) is considered with each element having a prescribed lifetime 
when subject to an applied stress (load). When elements fail, their load is redistributed to 
other elements of the set according to a prescribed rule of transfer. As a consequence 
of the load transfer, the lifetime of the receptors is actually reduced and the question is: 
how long does it take for the whole set to collapse? These fiber-bundle models are called 
dynamical, or time-dependent 0-§]], as opposed to their static counterparts, which have 
also been intensively studied The rule for redistributing the load of failed elements 

can be wide, but there are two limiting cases. In the first, the stress of the failed element 
is transferred equally to all surviving elements (ELS, for Equal Load Sharing). In the 
second, the load of the failed element is transferred to the nearest surviving element (s) 
(LLS, for Local Load Sharing). Hierarchically organized transfer (HLS) criteria are also of 
great interest Recently, these models have received much attention in the geophysical 



literature | ID| , because one would reasonably expect the emergence of universal scaling laws 
of the type observed in seismology |TT],|T2|. In this field, the bundle is a representation of a 
fault, and the individual elements or fibers represent asperities on the fault plane. 

In Reference [|L2 |, the ELS case is formulated in terms of a differential equation of the 
radioactive decay type. We have followed this perspective to devise a numerical probabilistic 
method to deal with any type of transfer rule. In Section |J, we explain in detail the 



differences between the usual approach and the new probabilistic approach. In Section |TJ, 
we compare results and present a brief discussion. 
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II. THE METHODS 



Suppose a set of N elements identified on the sites of a supporting lattice. This infor- 
mation is contained in a list {xi it } = {x}t , 1 < % < iV . This list is necessary, except in the 
ELS case, to know how to distribute the load of the failed elements. The broken elements 
are marked in {x} t - At t — 0, all the elements of the set are loaded with a reference value 
cr = 1. At any time, the total load acting on the surviving elements is constant, equal to 
N Q a . 

To each element, i, one assigns a random lifetime, ti^, under the unity of stress: 

ni = l-e"* ! - , (2.1) 

where rij are random numbers between and 1. This choice implies a logarithmic distribution 
of lifetimes. A more general distribution function for the failure time of a single element 
subjected to a known load history <r(t) is (see, for example @,[ll[) 



rii = 1 — exp 



(2.2) 



where \l/(x) is the shape function. The time integral in Eq. ( |2.2j ) introduces a hazard rate ^(cr) 
known generally as the breakdown rule in terms of the instantaneous load level. Experimental 
and theoretical work |2|,[|] favors a shape function \P(x) of the form 

^(x)=x /3 , (2.3) 

known as the Weibull shape function, with the particular choice (3 = 1 giving the exponential 
shape function. 

As for the breakdown rule, two special forms are widely used in the literature: the 
exponential breakdown rule, 

Ue = e ' ?(,T/,To) , (2.4) 

and the power-law breakdown rule, 
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v P = vo(-Y, (2.5) 

with 0, rj, p, u , o"o all positive constants, p is called the Weibull exponent because inserting 
Eq. ( |2.3| ) and Eq. ( |2.5| ) in Eq. (|2.2|) mimics the static Weibull distribution for the failure 
load of a single element. This parameter typically varies between 2 and 5. The exponential 
breakdown rule, Eq. ( |2.4j ) has a characteristic failure rate, whereas the power-law breakdown 



rule is scale free and can be regarded as a local approximation to the former. Following |12 
we will use Eq. Q2.5| ) as the individual breakdown rule in order to be able to compare 
the perfomance of the two approaches for load-transfer models. For further insights into 
the theoretical and experimental basis of the Weibull shape function see ||, and for the 
breakdown rules see References [f§] and ||TTf| . 



Without losing generality one can choose u = a = 1 in Eq. (|2.5| ), which means that a p 
is a measure of the failure rate (i.e, a unit failure rate under the unity of load is assumed). As 
vq is actually a frequency, z/ t is a dimensionless time variable and because of the particular 
choice Uq — 1, t will hereafter stand for non-dimensional time. 

If one substitutes Eq. ( |2.3| )and Eq. ( |2.5| ) in Eq. (|2.2j ) with the particular set of constants 
(3 = Uq = a = 1, we obtain 



71; 



1 - exp (- J l '° o p (T)dT^j , (2.6) 



which can be integrated for constant unit load cx(t) = <7q = 1 to give Eq. ( |2.1| ). 

When loads of failed elements are redistributed, the load acting on each element will no 
longer be the constant a but will depend on time o"j(t) > a = 1. Thus we introduce a 
reduced time to failure for each element, Tij, given by 

rTij 



ti,0 



^Ydt (2.7) 

In the case of independent elements, <7j(t) — ao — 1 and t^o = Tij. However, load transfer 
occurs, and hence the actual time to failure of element i, T t j is reduced to below £ i>0 . By 
imposing the fulfillment of Eq. ( |2.7| ), the successive order of breaking of the N elements, 



one after the other, is easily identified and the total time of collapse is the T;j of the longest 
lasting element. Thus, in this approach the randomness, that is the population of lifetimes, 
is fixed at t — (quenched disorder), and the breaking process is completely deterministic. 
Henceforth, we will refer to this approach as the usual one. 

In the new probabilistic approach presented here, the fluctuations are generated during 
the breaking process and hence it is an example of so-called annealed randomness. An 
interesting question is whether the two types of disorder, namely quenched and annealed 
disorder, in these models lead to different results, as has been observed for some critical 
phenomena M. 



In Ref. ||12|| , Newman et al formulated the ELS mode in terms of a differential equation 
of the radioactive decay type. Denoting the number of surviving elements as N s (t), its 
differential change is given by 

dN 

=T = - N °° P > ( 2 - 8 ) 
dt 

hence <j p represents the decay rate. But in the ELS mode o = hence 

N s (t) = N [p(T f -t)} 1 /" (2.9) 

and Tf = 1/p. In this setting, fluctuations do not exist and one simply obtains mean values 
for the failure rate. 

Following the perspective of the previous differential equation in which a group of ele- 
ments, supporting the same load a, fail at a rate a p or, in other words, have a mean-life 
l/(j p , one can devise a probabilistic method for any transfer rule. The scenario would be a 
set of iVo elements identified in {x} t , like a sample of radioactive nuclei fixed on a lattice, 
all having initially a decay rate aft = 1. As time passes, failures (disintegrations) occur and 
this does not merely imply the effective disappearance of the failed elements, but also the 
modification of the decay rate of other surviving elements. The modification comes from the 
redistribution of load as accorded in the rule of transfer (ELS, LLS, etc), and the assumption 
that the decay rate of any element is given by its a p value. As in this strategy of calculation 



one has to proceed at discrete time intervals, Sj, j = 1, 2, the information of loads in the 
set will be contained in a list denoted by {<Jij} = {&}j- The list is updated at each time 
step, together with {x}j. After j time steps, there will have appeared subsets of elements. 
Each subset is formed by all the surviving elements bearing the same load. We organise 
these subsets into sublists identified by the subindex I, and denote the corresponding load 
by Y\ and the number of elements belonging to the sublist / by N\. This information, which 
is obtained from {a}j, will be denoted by {Yi,Ni}j and updated simultaneously. At the 
beginning, as the load of all elements is 1, the sublists are 

^,0 = 1,^,0 = ^0 if/ = 1, 

Y lfi = 0,^,0 = if/ ^1. 

Now, it is clear that the simultaneous existence of several sublists in the sample, each 
with a different decay rate Yf, poses a difficulty for an accurate description of the decay 
process of the whole set The key point is the choice of the length of the time intervals, 
5j. To illustrate this problem, in Fig. [I] we have plotted the detailed evolution of breaking of 
a hierarchical set of iVo = 1024, coordination number equal to 2, and p = 4. In abscisas one 
represents time, from to Tf. In ordinates the spatial position of the N elements of the set 
is represented. At t = all elements are sound. As time evolves, breakings (represented by 
small crosses) are produced and therefore the number of failed elements, represented by the 
continuous line, grows. At t = Tf the number of failures is N . The height of the vertical 
spikes represents the load supported by an element at the time of failing. For short times, 
ruptures appear dispersed across the set and the rate of breaking is small. Progressively 
though, there appear cracks formed by the failures of neighboring elements, and this makes 
the continouos line adopt steeper slopes. Finally, the final breakdown occurs related to a 
big crack of a size similar to the whole system. This stage is also related to the high values 
of the spikes. The progressive acceleration of the breaking process is thus clear from this 
figure. 

Therefore the time interval used in the probabilistic approach must be variable with 
time. Otherwise, if one takes for 5 a reasonable value for the beginning, the final part of the 
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breaking will be badly described: in each 5, many elements will fail and the prediction of 
Tf will be very inaccurate. On the other hand, if one chooses a 5 typical of the final stages, 
it will be so short that although the calculation would be extremely good, at the beginning 
one would become bored awaiting the outcome of a breaking event. That is to say, these 
small intervals are not realistic for practical use. 

It is for this reason that when one tries to devise an efficient numerical method to 
accurately describe the time evolution of the system, the choice of the time interval must be 
adjusted to the characteristic scale at which individual elements break in the process. This 
characteristic time scale, 5, as mentioned before changes with time, 5j, and is implemented 
through the following definition 

5a = minimun of < — —t. — \ (2.10) 

where v is a constant > 1, independent of I and of j; we will call it the time resolution 
parameter. The length of 5j as defined in Eq. ( [2.10|) points, at each j, to a specific sublist 
whose / will be denoted as kj. Now we define a probability of failure for each sublist, 

PU = Y &j ■ ( 2 - n ) 

As Yf- is the failure rate for elements in sublist I, Yf^Sj is the expected number of 
casualties per element in sublist I, i.e, the probability of failure for sublist I. The product 
Ni j pij is maximum for the sublist kj, in this case N^ j p^j — which means that when 
the comparisons below (2.12) are perfomed, elements belonging to kj are the most likely to 



fail. In particular if v = 1, one element of the kth sublist is likely to fail. For the other 
sublists, the probability of an individual failure is lower than one. However, any element of 
any lj has a non-zero chance of failing in this probabilistic approach. We have called v the 
time resolution parameter because if it grows the time intervals 5j are smaller and therefore 
it is obvious that the process of failure is more finely resolved. Then the probability pij is 
compared, for each element belonging to the sublist with a random number, n, < n < 1. 

If pij > n the element fails. (2.12a) 



If pi j < n the element survives. 



(2.12b) 



The elements that fail in any of the sublists transfer their load according to the rule of 
transfer and the information contained in the list {x}j. In the case that no element fails, 
a new time interval 5j+i equal to 8j is added and the same probabilities, pij+i = pij, are 
compared with random numbers. This is repeated until at least one failure occurs which 
modifies {x}j, {o~}j, {Yi, Ni}j, and hence Sj. The total time to failure, Tf, is the sum of the 
Sj up to the disappearance of all the elements. 

In the ELS case, Sj can be explicitly written. After (j — 1) steps and assuming one failure 
per step, the number of surviving elements forming the unique sublist is Nj = Nq — (j — 1), 
and the individual load is = (N Q /Nj). Then 



Note that we have used v = 1 to be in accordance with the one failure per step assumption. 
In Eq. ( |2.13|) one observes that, in the first step, Si = l/N Q , and in the last step, Sn = l/^o- 
Now we proceed to sum up all the time intervals. In the continouos limit, we find 



For the ELS case, p = 2, Nq = 100, we plot in Fig. |^ the average of Tf after the number 
of simulations expressed in the abscisas, for various v. The horizontal lines comprise the 
extremes of the values obtained in 10 averages of 32000 simulations each by means of the 
usual method. One can observe, a) the actual Tf of this set is not - as predicted by the 
differential equation, this being a finite-size effect; and b) v = 4 is already sufficient in 
this method to reproduce the result of the usual approach. For the HLS case, p = 2 and 
coordination number of the Cayley tree, c = 2, we show in Fig. ^ the dispersion of Tf 




(2.13) 




(2.14) 



which tends to the correct result 1/p in the limit of large Nq. 



III. RESULTS AND DISCUSSION 



S 



emerging from the usual method (squares) and from the new method (circles) with v = 1. 
Note the slight shift rightwards of the center of the Gaussian, i.e, the values are longer. This 
is in agreement with what is seen in Fig. 0. A greater value of v would move the Gaussian to 
the left up to the coincidence. In Fig. |], the averaged rates of breaking of a set of Nq = 128 
are plotted under the HLS rule, c = 2, for two values of the Weibull exponent p = 4 and 
p = 6. We compare the habitual method and the new method for v = 1. In Table |, a set 
of values of Tf and their intrinsic width is shown for the HLS case, c = 2, by varying Nq 
and v. Data for iVo = 128, and Nq = 512 are averages over 32000, and 10000 realizations 
respectively. The errors quoted are one standard deviation of the mean. Perhaps the most 
abrupt rule of transfer that one can imagine is that of the local one- dimensional unilateral 
model IfR], where the load of failed elements is transferred to the nearest neighbor in the 
row going in one direction. This implies the almost immediate opening of big cracks and 
hence a great instability. The probabilistic approach for this model has been tested and 
again both methods coincide. 

Note that in the probabilistic method, there can be Sj in which no element fails, and 
others in which several elements do. In contrast to the usual method, here no disorder is fixed 
at the start: we begin with N elements, all with the same mean-life; the random successive 
failures are responsible for the fluctuations, i.e, this is an example of annealed disorder. 
For a small value of u, the results emerging from the probabilistic method are already 
indistinguishable from those deriving from the usual method. So, we have numerically 
proved the equivalence of the two approaches. If v is greater, the method demands more 
effort but the results reach a saturation point. Comparing the respective disadvantanges of 
computing: in the probabilistic approach it is necessary to deal with larger sets of random 
numbers while in the usual method the set of stored data is much bigger. 

We conclude by quoting Feynman who, in his original paper on path integral formalism 



I5|, writes "although it does not yield new results there is a pleasure in recognizing old 



things from a new point of view" . 
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FIGURES 

FIG. 1. Rate of breaking (continouos line) vs. time of a hierarchical set. Small crosses represent 
the position of local fractures. The height of the vertical spikes indicates the load of an element at 
the time of its failure. Read the text for details. 

FIG. 2. Tf results using probabilistic method for various v. The two horizontal lines show the 
prediction of the usual method. 

FIG. 3. Comparison of the two methods, for a hierarchical model. 

FIG. 4. Evolution of N s over time for both methods. 
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TABLES 



TABLE I. Comparison of the two methods (HLS, c=2)(see text for details). 





p = 2 




p = 4 




Method 


iVo=128 a 


iVo=512 b 


N =128 


iV =512 


Prob(> = 1) 


0.3890 


0.3614 


0.0992 


0.0750 


Prob(f = 4) 


0.3807 


0.3577 


0.0907 


0.0718 


Prob(^ = 10) 


0.3781 


0.3573 


0.0901 


0.0710 


Usual 


0.3776 


0.3565 


0.0888 


0.0705 



Simulations with iVo=128 elements are averages over 32000 realizations. Standard deviation of 
the mean value is ±2 units in the least significant digit. 

Simulations with iVo=512 elements are averages over 10000 realizations. Standard deviation of 
the mean value is ±1 unit in the least significant digit. 
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